Parametric  Reduced-Order  Models  for  Probabilistic 
Analysis  of  Unsteady  Aerodynamic  Applications 
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Methodology  is  presented  to  derive  reduced-order  models  for  large-scale  parametric 
applications  in  unsteady  aerodynamics.  The  specific  case  considered  in  this  paper  is  a 
computational  fluid  dynamic  (CFD)  model  with  parametric  dependence  that  arises  from 
geometric  shape  variations.  The  first  key  contribution  of  the  methodology  is  the  derivation 
of  a  linearized  model  that  permits  the  effects  of  geometry  variations  to  be  represented  with 
an  explicit  affine  function.  The  second  key  contribution  is  an  adaptive  sampling  method 
that  utilizes  an  optimization  formulation  to  derive  a  reduced  basis  that  spans  the  space  of 
geometric  input  parameters.  The  method  is  applied  to  derive  efficient  reduced-order  models 
for  probabilistic  analysis  of  the  effects  of  blade  geometry  variation  for  a  two-dimensional 
model  problem  governed  by  the  Euler  equations.  Reduced-order  models  that  achieve  three 
orders  of  magnitude  reduction  in  the  number  of  states  are  shown  to  accurately  reproduce 
CFD  Monte  Carlo  simulation  results  at  a  fraction  of  the  computational  cost. 


I.  Introduction 

Model  reduction  is  a  powerful  tool  that  permits  the  systematic  generation  of  cost-efficient  representations 
of  large-scale  systems  that,  for  example,  result  from  discretization  of  partial  differential  equations  (PDEs). 
We  present  an  approach  for  deriving  reduced  models  for  probabilistic  analysis  in  large-scale  unsteady  aero¬ 
dynamic  applications.  The  key  challenges  that  must  be  addressed  in  this  setting  are  the  formulation  of  a 
computationally  efficient  representation  of  the  parametric  dependence  that  describes  the  uncertainty,  and 
the  derivation  of  reduced-order  models  that  capture  variation  over  a  parametric  input  space. 

Recent  years  have  seen  the  widespread  use  of  computational  fluid  dynamics  (CFD)  to  solve  problems 
arising  in  applications  of  interest  in  unsteady  aerodynamics.  In  many  cases,  CFD  formulations  lead  to  large- 
scale  systems  of  equations  that  are  computationally  expensive  to  solve.  In  many  unsteady  aerodynamic 
applications,  a  small  number  of  inputs  and  outputs  of  interest  can  be  identified,  and  computationally  effi¬ 
cient  reduced-order  models  can  be  obtained  that  preserve  the  desired  input-output  mapping.  For  example, 
the  proper  orthogonal  decomposition  (POD)  method  of  snapshots1  has  been  used  widely  throughout  CFD 
applications  such  as  aeroelasticity2-4  and  flow  control.5,6 

Quantifying  the  impact  of  variations  in  input  parameters  on  system  outputs  of  interest  is  critical  to  a 
number  of  applications,  such  as  shape  design  and  probabilistic  analyses.  For  example,  mistiming,  or  blade- 
to-blade  variation,  is  an  important  consideration  for  aeroelastic  analysis  of  bladecl  disks,  since  even  small 
variations  among  blades  can  have  a  large  impact  on  the  forced  response  and  consequently  the  high-cycle 
fatigue  properties  of  the  engine.  In  such  applications — where  the  physical  system  must  be  simulated  repeat¬ 
edly  for  different  inputs — the  availability  of  reduced  models  can  greatly  facilitate  the  design  and/or  analysis 
task.  However,  to  be  useful  in  such  a  setting,  the  reduced  model  must  provide  an  accurate  representation  of 
the  high-fidelity  CFD  model  over  a  wide  range  of  parameters. 

Most  reduction  techniques  for  large-scale  systems  employ  a  projection  framework  that  utilizes  a  reduced- 
space  basis.  Methods  to  compute  the  basis  in  the  large-scale  setting  include  approximate  balanced  trunca¬ 
tion,7-10  Krylov-subspace  methods,11-13  and  POD.1,14, 15  In  the  latter  two  cases,  the  quality  of  the  reduced- 
order  model  is  critically  dependent  on  the  information,  generated  from  sampled  solutions  of  the  large-scale 
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system,  that  is  used  to  create  the  reduced  basis.  In  general,  selecting  an  appropriate  set  of  samples  to 
generate  this  information  has  been  via  an  ad-hoc  process.  Empirical  knowledge  of  the  problem  at  hand  has 
been  used  to  sample  a  parameter  space  to  generate  a  POD  or  Krylov  basis  for  cases  where  the  number  of 
input  parameters  is  small,  for  example  optimal  control  applications5,6,16,1'  and  parameterized  design  of  in¬ 
terconnect  circuits,18  and  in  the  case  of  multiple  parameters  describing  inhomogeneous  boundary  conditions 
for  parabolic  PDEs.19 

For  reduction  of  large-scale  linear  time-invariant  systems  using  multipoint  rational  Krylov  approxima¬ 
tions,  Ref.  20  proposes  a  systematic  method  for  selecting  interpolation  points  based  on  a  rigorous  optimality 
criterion.  To  address  the  more  general  challenge  of  sampling  a  high-dimensional  parameter  space  to  build  a 
reduced  basis,  the  greedy  algorithm  was  introduced  in  Refs.  21-24.  The  key  premise  of  the  greedy  algorithm 
is  to  adaptively  choose  samples  by  finding  the  location  of  maximum  reduced  model  error  estimate,  over  a 
pre-determined  discrete  set  of  parameters.  The  greedy  algorithm  was  applied  to  find  reduced  models  for 
the  parameterized  steady  incompressible  Navier-Stokes  equations.22  It  was  also  combined  with  a  posteri¬ 
ori  error  estimators  for  parameterized  parabolic  PDEs,  and  applied  to  several  optimal  control  and  inverse 
problems.23,24 

Here,  we  formulate  the  task  of  determining  appropriate  sample  locations  as  a  greedy  optimization  prob¬ 
lem,  which  is  solved  using  an  efficient  adaptive  algorithm.  The  optimization  formulation  treats  the  parameter 
space  as  continuous;  that  is,  we  do  not  require  the  a  priori  selection  of  a  discrete  parameter  set.  Further,  our 
selection  criteria  uses  the  true  computed  error  between  full-order  and  reduced-order  outputs;  thus,  our  ap¬ 
proach  is  applicable  in  cases  for  which  error  estimators  are  unavailable.  Unlike  other  sampling  methods,  the 
optimization-based  approach  scales  well  to  systems  with  a  large  number  of  parameters.  To  further  address 
the  challenge  of  achieving  a  computationally  efficient  representation  of  the  dependence  of  the  CFD  model 
on  geometric  parameters,  we  propose  a  linearization  strategy  that  yields  an  affine  parametric  dependence. 

This  article  is  organized  as  follows.  Section  II  describes  the  discontinuous  Galerkin  (DG)  CFD  model  used 
in  this  work  and  formulates  a  linearized  model  for  capturing  the  effects  of  geometry  variations  on  unsteady 
aerodynamic  response.  Section  III  presents  an  overview  of  projection-based  model  reduction  and  then 
describes  the  proposed  optimization-based  approach  to  determine  the  reduced  basis.  Section  IV  demonstrates 
the  methodology  through  an  example  that  considers  the  effects  of  variations  in  blade  geometry  on  the  forced 
response  of  a  subsonic  compressor  blade  row.  Finally,  Section  V  presents  conclusions. 

II.  Linearized  CFD  Model  with  Geometric  Variability 

In  this  section,  we  present  an  overview  of  the  DG  CFD  model  that  is  used  in  this  work.  We  then  derive 
a  linearized  unsteady  model  that  incorporates  geometric  variability  in  a  computationally  efficient  way. 

A.  CFD  model 

Recent  developments  in  the  field  of  CFD  have  led  to  the  use  of  higher-order  finite  element  discretizations  for 
flow  modeling.  These  schemes  have  advantages  over  traditional  finite- volume  methods  by  introducing  higher- 
order  accuracy  compactly  within  grid  elements  and  thus  providing  a  significant  decrease  in  the  computational 
cost  to  obtain  reliably  accurate  solutions.  A  DG  formulation  is  used  in  this  work.  A  linearized  unsteady 
flow  solver  was  developed  and  implemented  as  part  of  a  larger  CFD  development  effort  that  includes  an 
adaptive  meshing  utility,  a  multigrid  solution  algorithm,  gradient-based  optimization  capability,  and  high- 
order  visualization.25  Here,  we  briefly  review  the  DG  discretization  and  solution  method  for  the  two- 
dimensional  Euler  equations,  which  are  described  in  more  detail  in  Ref.  26  and  27. 

1.  Nonlinear  CFD  model 

For  two-dimensional  flows,  the  Euler  equations  are  given  by 

^+V-^(w)  =  0,  (1) 

where  w  =  [p,  pu,  pv,  pE]T  is  the  conservative  state  vector,  and  T  =  (Fa:,Fy)  is  the  inviscid  Euler  flux, 
where  Fx  =  [ pu ,  pu2  +  P,  puv,  puIl\T  and  Fv  =  [pv,  puv,  pv2  +  P,  pvH]T .  Here,  p  is  the  density,  u  and 
v  are  respectively  the  x —  and  y— component  of  velocity,  E  is  the  energy,  P  is  the  pressure,  and  H  =  E  +  P/p 
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is  the  total  enthalpy.  The  equation  of  state  is 

P=(7-l)  pE-^p(u2  +  v2)  ,  (2) 

where  7  is  the  ratio  of  specific  heats. 

As  in  the  continuous  finite  element  method,  the  first  step  in  the  DG  method  is  to  discretize  the  domain 
under  consideration,  Q,  into  elements  f le.  Next,  a  space  of  polynomials  of  degree  at  most  p,  is 

defined  on  each  element,  where  h  denotes  a  representative  element  size  for  the  discretization  (e.g.  the  size 
of  the  smallest  element).  On  each  element  f le,  the  approximate  solution  w/,  can  be  found  by  enforcing  the 
nonlinear  conservation  law  (1)  locally,  for  all  test  functions 

[  -  I  Vv[T(w h)dne+  [  (v+)TW(w+,w^,n)ds 

Jne  01  Jon*  Jan^d  n 

+  [  (vh)T^bd(w+,w-,n)ds  =  0,  (3) 

Janendn 

where  dfl  and  d are  the  boundaries  of  the  entire  domain  O  and  the  element  f le,  respectively,  and  n 
denotes  the  outward-pointing  normal  on  the  boundaries  of  the  element.  The  terms  7Y(w0,  w|0  n)  and 
Wbd(w0,w0,n)  are  numerical  flux  functions  for  interior  and  boundary  edges,  respectively,  where  ()+  and 
()-  denote  values  taken  from  the  interior  and  exterior  of  the  element.  The  interior  flux  function  is  computed 
using  the  Roe-averaged  flux  function28  and  contributes  over  element  boundaries  that  do  not  belong  to  the 
domain  boundary,  denoted  by  d£le\dQ.  The  fluxes  on  the  common  boundaries  of  dCle  and  df 1,  denoted  by 
dfle  fl  dfi,  are  computed  using  the  inner  state  and  boundary  condition  data. 

The  final  form  of  the  DG  discretization  is  constructed  by  selecting  a  basis  for  The  approximate 

solution  w;,  on  each  element  is  assumed  to  be  a  linear  combination  of  the  finite  element  basis  functions  if>j, 

nb 

w  h(t,x,y)  =  (4) 

j= 1 

where  w.;  (t)  gives  the  modal  content  of  ifij  on  element  f le,  and  rii,  is  the  number  of  basis  functions  required 
to  describe  (e.g.  rib  =  1  for  p  =  0  and  rib  =  3  for  p  =  1).  The  complete  set  of  unknown  quantities 

for  the  DG  formulation  thus  comprises  the  values  of  w j(t),  j  =  1, . . .  ,nb  for  every  element  in  the  spatial 
domain.  These  quantities  are  contained  in  the  vector  w(t)  £  IR",  where  n  is  the  total  number  of  unknowns, 
which  depends  both  on  the  number  of  elements  in  the  discretization  and  on  the  polynomial  order  p. 

This  spatial  discretization  together  with  the  application  of  appropriate  boundary  conditions  leads  us  to 
the  following  set  of  nonlinear  ordinary  differential  equations  (see  Ref.  25  and  27  for  more  detail), 

E^  +  R(w,U)=°,  (5) 

where  E  £  ]Rraxn  is  the  mass  matrix,  and  R  is  the  residual  vector  representing  the  final  three  terms  of  (3) 
plus  the  effects  of  boundary  conditions.  The  vector  u(t)  £  JR"'  contains  m  external  forcing  inputs  that  are 
applied  through  boundary  conditions,  such  as  prescribed  motion  of  the  domain  boundary  or  incoming  flow 
disturbances. 

For  steady-state  flows,  pseudo  time-stepping  is  used  to  improve  the  initial  transient  behavior  of  the 
solver  and  the  nonlinear  system  (5)  is  solved  using  a  p-multigrid  scheme  with  a  line  Jacobi  smoother.25,27 
For  unsteady  computations,  a  second-order  backward  Euler  temporal  discretization  is  applied  to  (5).  The 
resulting  nonlinear  equations  are  solved  using  a  Newton  solver.  The  motion  of  grid  points  on  the  domain 
boundary  is  prescribed  according  to  the  corresponding  external  input  (e.g.  prescribed  motion  of  a  blade). 
The  resulting  motion  of  internal  grid  points  is  computed  using  a  Jacobi  smoothing  formulation. 

For  many  practical  applications,  we  are  concerned  with  the  prediction  of  a  set  of  l  output  quantities  of 
interest.  We  define  these  output  quantities  of  interest  to  be  contained  in  the  vector  y  €  10  and  defined  by 
the  nonlinear  function  Y, 

y(t)=Y(w(t)).  (6) 
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2.  Linearized  CFD  model 


In  many  cases  of  interest,  the  unsteady  flow  solution  can  be  assumed  to  be  a  small  perturbation  from  steady- 
state  conditions.  This  allows  the  unsteady  governing  equations  to  be  linearized  around  the  steady-state  flow, 
which  reduces  the  computational  cost  of  solution  considerably.  The  linearization  of  equations  (5,  6)  can  be 
written  in  standard  state-space  form 


Ex  =  Ax  +  Bu,  (7) 

y  =  Cx,  (8) 


where  x  £  3Rra  is  the  state  vector  containing  the  n  perturbations  in  flow  unknowns  from  the  steady-state 
solution  wss,  that  is  w(t)  =  wss  +  x(i).  The  matrices  A  €  ]Rnxn,  B  £  JRnxp,  and  C  £  Mgxn  in  (7)  and  (8) 
have  constant  coefficients  evaluated  at  steady-state  conditions  and  arise  from  the  linearization  of  (5)  and  (6) 
as  follows, 


A  = 


<9R 

dw 


B  = 


dR 

du 


C  = 


dY 

dw 


(9) 


B.  Linearized  unsteady  model  with  geometric  variability 

We  next  consider  incorporating  the  effects  of  geometry  variability  into  the  linearized  unsteady  CFD  model. 
Following  Ref.  29,  a  general  geometry,  g,  can  be  expressed  as 

ns 

g  =  9n  + 9 +  ^2^iZiVi,  (10) 

i=l 

where  g„  is  the  nominal  geometry,  g  is  the  average  geometric  variation,  ty  are  geometric  mode  shapes,  and 
ns  is  the  number  of  mode  shapes  used  to  represent  the  variation  in  geometry.  The  geometric  mode  shapes 
could  be  computed,  for  example,  by  performing  the  principal  component  analysis  (PCA)  on  a  manufacturing 
sample  of  system  geometries.  In  that  case,  the  parameters  zt  in  (10)  are  random  numbers  normally  distributed 
with  zero  mean  and  unity  variance,  Zi  £  N(0, 1),  and  <7,;  is  the  standard  deviation  of  the  geometric  data 
attributable  to  the  ith  mode;  thus  the  product  aiZi  is  the  amount  by  which  the  mode  vt  contributes  to  the 
geometry  g.  A  detailed  description  of  the  methodology  underlying  this  geometric  model  can  be  found  in 
Ref.  29. 

Using  the  model  (10),  a  general  geometry  g( z)  is  specified  by  the  parameter  vector  z  =  [zi,  Z2,  ■  ■  ■ ,  zn JT, 
which  describes  the  geometry  variability  in  terms  of  the  geometry  modes.  The  linearized  CFD  system 
corresponding  to  geometry  g{ z)  is  given  by 

E(wss(c/(z)),  g(z))x  =  A(wss(ff(z)),flr(z))x  +  B(wss(flr(z)),5(z))u,  (11) 

y  =  C(wss(g(z)),g(z))x,  (12) 

where  the  CFD  system  matrices  E,  A,B  and  C  are  in  general  both  a  nonlinear  function  of  the  geometry, 
g( z),  and  of  the  steady-state  solution,  wss(c/(z)),  which  is  itself  also  a  function  of  the  geometry.  To  solve  the 
CFD  system  (11),  (12),  for  each  geometry  g  we  must  firstly  compute  the  steady-state  solution,  wss(g(z)), 
secondly  evaluate  the  linearized  matrices  E,A,B  and  C,  and  thirdly  solve  the  resulting  large-scale  linear 
system.  This  is  a  computationally  prohibitive  proposition  for  applications  such  as  probabilistic  analysis, 
where  thousands  of  geometry  perturbations  may  be  analyzed  over  many  random  samples  z. 

For  convenience  of  notation,  we  write  the  dependence  of  the  CFD  matrices  on  the  parameter  z  as 
E(wss(g(z)),g(z))  =  E(z),  A(wss(g(z)),  g(z))  =  A(z),  B{waa(g(z)),g(z))  =  B(z),  and  C(wss(g(z)),  g(z))  = 
C(z).  We  use  the  expansion  given  by  equation  (10),  which  represents  a  general  geometry  as  a  perturbation 
about  the  average  geometry  go  =  gn  +  g,  to  derive  an  approximate  model  for  representing  the  effects  of 
geometry  variations.  Instead  of  computing  the  linearized  CFD  matrices  exactly  for  any  random  variability 
z,  we  choose  to  linearize  the  relationships  E(z),  A(z),  B(z),  and  C(z).  We  define  the  linearized  unsteady 
CFD  model  for  the  average  geometry  go  =  gn  +g  by  the  matrices  Eo,  Ao,  Bo,  and  Co,  with  corresponding 
solution  xq.  That  is,  for  z  =  0we  have 


E0x0  =  A0x0  +  B0u, 


y0  =  C0x0. 


(13) 

(14) 
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Using  a  Taylor  series  expansion  about  z  =  0  for  the  matrix  A(z)  gives 


A(z)  =  An 


<9A 

dz^ 


z  1 


dA 


z— 0 


z— 0 


(15) 


where  the  matrix  partial  derivatives  denote  componentwise  derivatives,  which  can  be  evaluated  through 
application  of  the  chain  rule.  These  derivatives  are  evaluated  at  average  geometry  conditions,  z  =  0.  The 
matrices  E(z),B(z)  and  C(z)  can  be  expanded  using  formulae  analogous  to  (15). 

If  the  geometric  variability  (given  by  the  product  (JiZi)  is  sufficiently  small,  the  constant  and  linear  terms 
in  the  Taylor  expansion  (15)  are  sufficient  to  approximate  the  linearized  matrices  A(z)  accurately,  that  is, 


A(z) 


Aq  + 


dA 

dzi 


Z\  +  . . .  + 

z— 0 


dA 

dzns 


Z=0 


For  i  =  1,2,...,  na,  we  define 


<9E 

dA 

-  <9B 

dc 

Ei  =  hU 

OZi 

5 

z=0 

z— 0 

dzi 

? 

z—0 

(16) 


(17) 


where  the  matrices  Ei,  A j,  B,,  and  C;  can  be  computed,  for  example,  using  a  finite  difference  approximation 
of  the  respective  derivatives.  The  approximate  linearized  CFD  model  for  any  geometric  variability  z  is  then 
given  by 


E0  +  ^  E jZj  ]  x  =  (  A0  +  ^  AjZj  )  x  +  (  B0  +  ^  B jZj  )  u, 


i- 1 


E(z) 


i—1 


i= 1 


y  = 


A(z) 

f  ns  \ 

Co  +  'y '  c iZi  I  x. 

V,  i= 1  ) 


B(z) 


(18) 


(19) 


C(z) 


It  should  be  noted  here  that  a  number  of  large-scale  steady-state  CFD  solves  are  required  in  order  to 
determine  the  matrices  Ao,  Bo,  Co,  Eo,  A*,  B;,  Cj  and  Ej.  For  example,  if  central  difference  approximations 
to  the  matrix  derivatives  are  used,  a  total  of  2  X)"=i  +1  large-scale  steady-state  CFD  solves  is  required.  This 
is  a  one-time  offline  cost;  once  the  matrices  are  computed,  the  approximate  linearized  system  (18),  (19)  can 
be  readily  evaluated  for  an  arbitrary  geometry  g{ z)  without  running  the  CFD  steady  solver. 

It  should  also  be  noted  that  the  model  (18),  (19)  is  valid  only  for  small  variations  from  the  average 
geometry.  Larger  variations  will  incur  larger  errors,  due  to  the  neglect  of  the  higher-order  terms  in  the 
Taylor  series  expansion.  Even  with  this  restriction,  the  model  is  useful  for  many  applications  where  small 
geometric  variations  are  of  interest;  however,  the  approximate  linearized  model  is  still  of  high  dimension,  and 
thus  is  computationally  too  expensive  for  applications  such  as  probabilistic  analysis  in  which  one  needs  to 
determine  the  unsteady  aerodynamic  response  for  many  random  geometries.  In  the  next  section  we  propose 
a  model  reduction  method  that  enables  us  to  further  reduce  the  cost  of  solving  the  approximate  linearized 
system.  The  key  challenge  that  must  be  addressed  is  developing  a  reduced-order  model  that  is  accurate  over 
both  time  and  the  geometric  parameter  space,  described  here  by  the  vector  z. 


III.  Model  Reduction  Methodology 

A.  General  projection  framework 

Most  large-scale  model  reduction  frameworks  are  based  on  a  projection  approach,  which  can  be  described  in 
general  terms  as  follows.  Although  the  projection  framework  is  described  here  for  a  general  system  that  is 
linear  in  the  state  variables,  it  can  also  be  applied  to  nonlinear  systems.  Consider  the  general  parameterized 
dynamical  system 


(20) 
(21) 
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with  initial  condition 


x(0)  =  x°,  (22) 

where  x(z ,t)  £  R"  is  the  state  vector,  u(t)  £  Rm  contains  the  m  forcing  inputs  to  the  system,  y(z ,x,f)  £ 
R ( .  contains  the  l  outputs  of  interest,  and  x°  is  the  specified  initial  state.  The  matrices  E  £  Rraxn, 
A  £  Rrax",  B  <E  R"xm,  and  C  £  Rix™  in  (20)  and  (21)  may  depend  (possibly  nonlinearly)  on  a  set  of  ns 
parameters  contained  in  the  vector  z  €  R"s.  General  dynamical  systems  of  the  form  (20)-(22)  often  arise 
from  discretization  of  PDEs.  In  that  case,  the  dimension  of  the  system,  n,  is  large,  and  the  parameters  z, 
could  describe,  for  example,  changes  in  the  domain  shape,  boundary  conditions,  or  PDE  coefficients.  The 
CFD  model  (18),  (19)  derived  in  the  previous  section  is  one  example  of  a  system  of  the  form  (20),  (21);  in 
that  case  the  parameters  Zj  describe  geometric  shape  variations. 

A  reduced-order  model  of  (20)-(22)  can  be  derived  by  assuming  that  the  state  x(z ,t)  is  represented  as  a 
linear  combination  of  nr  basis  vectors, 

x  =  €>xr,  (23) 

where  x(z ,t)  is  the  reduced  model  approximation  of  the  state  x(z ,t)  and  nr  <C  n.  The  projection  matrix 
<1?  £  Rnxrar  contains  as  columns  the  basis  vectors  fa,  i.e.,  <l>  =  [fa  fa  •  •  •  </>„,.],  and  the  vector  xr(z,  t)  £  E"r 
contains  the  corresponding  modal  amplitudes.  Using  the  representation  (23)  together  with  a  Petrov-Galerkin 
projection  of  the  system  (20)— (22)  onto  the  space  spanned  by  the  left  basis  \F  £  R”xriT'  yields  the  reduced- 
order  model  with  state  xr(z ,t)  and  output  yr(z,x,f), 


Er(z)xr  =  Ar(z)xr  +  Br(z)u, 

(24) 

yr  =  Cr(z)xr, 

(25) 

1  o 

II 

£ 

o 

(26) 

where  Er(z)  =  ^TE(z)^,  Ar(z)  =  tfTA(z)$,  Br(z)  =  ’®fTB(z),  and  Cr(z)  =  C(z)$. 

Projection-based  model  reduction  techniques  seek  to  find  the  bases  and  such  that  the  reduced 
system  (24) — (26)  provides  an  accurate  representation  of  the  large-scale  system  (20)-(22)  over  the  desired 
range  of  inputs  u (t)  and  parameters  z.  Here,  we  consider  Galerkin  projections,  that  is  4?  =  4/ ,  although  our 
methodology  holds  in  the  general  case. 

B.  Reduced  basis  for  parametric  input  dependence 

Using  the  general  projection  framework,  our  model  reduction  task  becomes  one  of  determining  an  appropriate 
reduced  basis  that  spans  both  the  parametric  input  space  z  and  the  space  of  unsteady  inputs  u(i).  In  the  case 
of  a  linear  time-invariant  system,  that  is,  a  system  of  the  form  (20) — (22)  with  no  dependence  on  parameters 
z,  a  number  of  model  reduction  techniques  can  be  used,  such  as  Krylov-based  methods  and  POD.  To  extend 
these  techniques  to  the  general  case  where  the  system  matrices  depend  on  the  parameters  z,  we  require  a 
systematic  method  of  sampling  the  parametric  input  space. 

In  the  case  of  the  POD,  the  reduced  basis  is  formed  as  the  span  of  a  set  instantaneous  state  solutions, 
commonly  referred  to  as  snapshots.  These  snapshots  are  computed  by  solving  the  system  (20)-(22)  for 
selected  values  of  the  parameters  z  and  selected  forcing  inputs  u(t).  The  quality  of  the  resulting  reduced- 
order  model  is  very  dependent  on  the  choice  of  parameters  and  inputs  over  which  snapshots  are  computed. 
Two  issues  arise  in  selecting  an  appropriate  sample  set.  First,  choosing  where  and  how  many  samples  to 
generate  is,  in  general,  an  ad-hoc  process.  One  can  use  knowledge  of  the  application  at  hand  to  determine 
representative  inputs;  however,  there  exist  no  guarantees  on  the  quality  of  the  resulting  reduced-order  model. 
Second,  in  the  case  that  the  parametric  input  space  is  of  high  dimension,  the  number  of  high-fidelity  system 
solves  required  to  generate  the  snapshots  in  an  ad-hoc  manner  can  become  prohibitively  large.  Using  standard 
sampling  methods,  a  problem  with  just  a  few  parameters  can  require  a  large  number  of  samples  to  adequately 
cover  the  space,  due  to  the  combinatorial  explosion  of  the  number  of  possible  parameter  combinations. 

Here,  we  use  the  greedy  algorithm2 1-24  to  adaptively  select  snapshots,  by  finding  the  location  in  parameter 
space  where  the  error  between  the  full-order  and  reduced-order  models  is  maximal,  updating  the  basis  with 
information  gathered  from  this  sample  location,  forming  a  new  reduced  model,  and  repeating  the  process. 
We  formulate  the  greedy  approach  as  an  optimization  problem  that  targets  the  error  in  reduced  model  output 
prediction,  which  is  defined  by  introducing  as  constraints  the  systems  of  equations  representing  the  full  and 
reduced  models.  The  optimization  formulation  treats  the  parameter  space  as  continuous;  that  is,  we  do  not 
require  the  a  priori  selection  of  a  discrete  parameter  set.  Further,  since  the  optimization  problem  uses  the 
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true  computed  error  between  full  and  reduced  outputs,  our  approach  is  applicable  in  cases  for  which  error 
estimators  are  unavailable. 

On  each  cycle  of  the  greedy  algorithm,  the  key  step  is  to  determine  the  location  in  parameter  space  where 
the  error  in  the  reduced  model  is  maximal.  We  define  the  cost  functional 

g(z  ,x,xr)  =  ^y  ||y(z,x,£)  -  yr(z,xr,f)|||  dt=^J  ||Cx  -  Crxr||  %  dt,  (27) 

which  describes  the  error  between  the  full  and  reduced  models  over  the  parameter  space  z,  integrated  over 
some  time  horizon  of  interest  tf.  Given  a  current  basis  <1?,  we  find  the  location  in  parameter  space  of 
maximum  error  by  solving  the  optimization  problem 


subject  to 


.  1 
max  y  =  — 

x,xr,z  2 

[  Cx  —  Crxr  |  dt 

Jo 

(28) 

E(z)x 

=  A(z)x  +  B(z)u, 

(29) 

x° 

=  x(0), 

(30) 

Er  (z)xr 

=  Ar(z)xr  +  Br(z)u, 

(31) 

x° 

II 

(32) 

^ min  — 

Z  ^  Z maxi 

(33) 

where  z min  and  zmax  are  respectively  lower  and  upper  bounds  on  the  parameter  vector  z.  We  denote  the 
parameter  vector  that  solves  the  maximization  problem  (28) — (33)  by  z*.  Next,  we  compute  the  solution 
x(z *,f)  of  the  full  system  at  the  worst-case  parameter  value  z*.  This  solution  information  is  added  to  the 
basis  3>,  for  example  using  the  POD.  (Note  that  once  the  sample  location  has  been  found,  other  model 
reduction  methods  could  also  be  employed.)  The  procedure  is  then  repeated  by  solving  the  optimization 
problem  (28)-(33)  with  the  updated  basis  3>.  Thus,  we  are  using  a  systematic,  adaptive  error  metric  based 
on  the  ability  of  the  reduced-order  model  to  capture  the  outputs  of  interest  in  order  to  choose  the  snapshot 
locations.  This  model  reduction  approach  is  summarized  in  the  following  algorithm. 

Algorithm  1  Adaptive  Sampling  Procedure 

1.  Given  a  reduced  basis  <l>.  solve  the  optimization  problem  (28) -(33)  to  find  the  location  in  parameter 
space  at  which  the  error  is  maximized ,  i.e.  find  z*  =  argmaxt/(z). 

2.  If  Q{ z*)  <  e,  where  e  is  the  desired  level  of  accuracy,  then  terminate  the  algorithm.  If  not,  go  to  the 
next  step. 

3.  With  z  =  z* ,  solve  the  full  system  (20)-(22)  to  compute  the  state  solutions  x(z *,£),  t  =  (0 ,£/).  Use 
the  span  of  these  state  solutions  to  update  the  basis  3>.  Go  to  step  1. 

In  Step  3  of  Algorithm  1,  the  basis  can  be  updated  using  many  of  the  existing  model  reduction  methods. 
For  example,  the  POD  could  be  used  to  compute  the  span  of  the  updated  snapshot  set,  which  would  comprise 
the  existing  basis  vectors  and  the  new  state  solutions  x(z *,£).  As  an  alternative  approach,  one  could  also 
solve  an  (inner)  optimization  problem  to  find  the  basis  that  minimizes  the  output  error  at  the  sample  points.30 
Algorithm  1  is  initialized  by  choosing  the  initial  basis  as  the  empty  set,  =  0;  thus  the  reduced  model  is 
initially  a  zero-order  approximation  of  the  full  model. 

The  optimization  problem  (28)-(33)  that  must  be  solved  in  each  adaptive  cycle  (i.e.  Step  1  of  Algorithm  1) 
is  large  scale;  in  particular,  note  that  the  large-scale  system  equations  appear  as  constraints  in  (29).  The 
determination  of  each  sample  point  z*  via  solution  of  this  optimization  problem  therefore  requires  some 
number  of  solves  of  the  system  (29),  which  for  the  large-scale  problems  of  interest  here  (n  >  105)  is  the 
dominant  computational  cost.  It  is  therefore  critical  to  use  an  efficient  optimization  method;  that  is,  one 
that  exploits  the  structure  of  the  problem  to  offer  rapid  convergence.  We  employ  recent  advances  in  scalable 
algorithms  for  large-scale  optimization  of  systems  governed  by  PDEs,  which  have  permitted  solution  of 
problems  with  millions  of  state  and  optimization  variables.31,32 
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In  order  to  solve  the  constrained  optimization  problem  (28)-(33),  we  choose  to  solve  an  equivalent  bound- 
constrained  optimization  problem  in  the  z  variables  by  eliminating  the  state  variables  x  and  xr.  That  is,  we 
replace  minx  x  z  <?(x,  xr,  z)  with  minz  t?(x(z),  xr(z),  z),  where  the  dependence  of  x  and  xr  on  z  is  implicit 
through  the  full  and  reduced  state  equations  (29)-(32).  To  deal  with  the  bound  constraints,  we  use  the 
Coleman-Li  approach  (see,  for  example  Ref.  33).  This  enables  us  to  use  the  subspace  trust-region  interior- 
reflective  Newton  framework,  proposed  in  Ref.  33,  to  solve  the  resulting  bound-constrained  optimization 
problem  efficiently.  That  is,  we  use  the  conjugate  gradient  (CG)  method  to  determine  the  subspace  over 
which  the  linear  system  of  equations  arising  at  each  Newton  step  is  solved,  and  globalize  by  a  trust  region 
scheme  (see,  for  example,  Ref.  34).  This  method  combines  the  rapid  locally-quadratic  convergence  rate 
properties  of  Newton’s  method,  the  effectiveness  of  trust  region  globalization  for  treating  ill-conditioned 
problems,  and  the  Eisenstat- Walker  idea  of  preventing  oversolving.35 

The  gradient  of  Q  with  respect  to  z,  as  required  by  Newton’s  method,  can  be  computed  efficiently  by 
an  adjoint  method.  The  Hessian-vector  product  as  required  by  CG  is  computed  on-the-fly;  because  it  is  a 
directional  derivative  of  the  gradient  its  computation  similarly  involves  solution  of  state-like  and  adjoint-like 
equations.  Therefore,  the  optimization  algorithm  requires  solution  of  a  pair  of  state  and  adjoint  systems  at 
each  CG  iteration. 

Since  the  system  dependence  on  the  parameter  z  is  nonlinear,  in  the  general  case  the  optimization  problem 
(28)— (33)  is  non-convex.  In  particular,  as  the  greedy  algorithm  progresses  we  expect  the  cost  functional  to 
become  increasingly  multi-modal,  since  the  error  function  will  be  close  to  zero  (below  the  tolerance  e)  at 
each  of  the  previous  parameter  sample  locations.  It  should  be  noted  that,  while  finding  the  global  maximum 
is  obviously  preferred,  convergence  to  a  local  maximum  is  not  necessarily  an  adverse  result.  Solving  the 
greedy  optimization  problem  is  a  heuristic  to  systematically  find  “good”  sample  points;  at  a  local  maximum 
the  error  is  (locally)  large.  To  avoid  convergence  to  a  local  maximum  close  to  a  previous  sample  location, 
and  thus  explore  the  parameter  space  more  widely,  a  random  initialization  of  the  optimization  variables  z 
is  used  for  each  cycle  of  Algorithm  1.  An  initial  guess  is  accepted  only  if  it  is  sufficiently  far  from  previous 
sample  locations,  measured  using  a  tolerance  that  is  set  relative  to  the  parameter  ranges.  The  stopping 
criterion  applied  in  Step  2  of  Algorithm  1  monitors  Q( z*),  the  reduced  model  error  associated  with  the 
optimal  solution  z*.  It  is  important  to  note  that  if  z*)  falls  below  the  desired  error  level,  this  guarantees 
only  that  the  local  error  between  full  and  reduced  model  is  sufficiently  small.  Due  to  the  non-convexity  of 
the  optimization  problem,  it  is  possible  that  larger  errors  may  exist  elsewhere  in  the  parameter  space. 

C.  Reduced-order  linearized  aerodynamic  model  with  geometric  variability 

Combining  the  linearized  unsteady  model  with  geometric  variability  from  Section  II  together  with  the  reduced 
basis  model  reduction  methodology  based  on  adaptive  sampling,  we  now  have  a  method  to  create  efficient 
reduced-order  models  that  capture  the  effects  of  small  geometric  variations. 

Using  the  projection  framework,  and  a  basis  <1?  computed  using  Algorithm  1,  the  reduced-order  model  of 
(18),  (19)  is 


Ero  +  ^  ^  E ViZi  j  xr  —  I  Aro  H-  ^  ^  AriZi  J  xr  H-  f  Bro  H-  ^  ^  B riz%  j  u, 


(34) 


i= 1 


Er(z 


i= 1 


i= 1 


Ar(z 


Br(z 


yr  =  C 


r*o  1  ^  'JGjriZi  J  Xr, 

i=  1  J 


(35) 


Cr(z 


where  the  reduced-order  matrices  are  given  by 

Ero  =  $tE0$,  Aro  =  $TAn$,  Bro  =  $TB0,  Cr„  =  C0$, 

Er.  =  4>tE,<I>,  Ar.  =  <3>7  A,<f>,  Br.  =  $tB},  Cr.  =  C,4>,  i  =  1, 2, . . . , ns. 


(36) 

(37) 


The  key  enabling  feature  of  the  adaptive  sampling  approach  is  that  it  allows  the  basis  <1?  to  be  computed  in 
an  efficient  systematic  manner,  even  when  the  dimension  of  the  parameter  space  is  large.  The  methodology 
also  gives  us  a  means  to  monitor  the  (local)  error  between  reduced-order  and  full-order  outputs.  The  key 
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advantage  of  the  linearized  variability  model  is  that  it  leads  to  an  affine  parameter  dependence;  thus  the 
reduced-order  matrices  in  (36)  and  (37)  can  be  evaluated  offline,  and  the  online  cost  of  solving  the  reduced- 
order  model  (34),  (35)  does  not  depend  on  the  large-scale  state  dimension  n. 

IV.  Probabilistic  Analysis  Application 

The  model  reduction  methodology  is  applied  to  probabilistic  analysis  of  a  subsonic  rotor  blade  that 
moves  in  unsteady  rigid  motion.  The  analysis  seeks  to  quantify  the  effects  on  the  blade  forced  response  of 
small  variations  in  blade  geometry.  Mistiming,  or  blade-to-blade  variation,  is  an  important  consideration 
for  aeroelastic  analysis  of  bladed  disks,  since  even  small  variations  among  blades  can  have  a  large  impact  on 
the  forced  response  and  consequently  the  high-cycle  fatigue  properties  of  the  engine.  The  effects  of  blade 
structural  mistiming  (variations  in  mass  and  stiffness  properties)  have  been  extensively  studied,  see  for 
example  Refs.  36-39;  however,  due  to  the  prohibitively  high  computational  cost  of  performing  probabilistic 
analysis  with  a  CFD  model,  the  aerodynamic  effects  due  to  variations  in  geometry  are  less  understood. 

Geometric  mistiming  effects  have  been  incorporated  into  structural  responses  of  bladed  disks  using  a 
mode-acceleration  method  to  convert  the  effect  of  geometric  mistiming  to  that  of  external  forces  of  the 
tuned  disks.40  Truncated  sets  of  tuned  system  modes  compensated  by  static  modes — generated  by  external 
forces  that  were  constructed  from  mistiming — were  then  used  to  obtain  efficient  and  accurate  structural 
reduced  models.  Several  studies  have  also  found  that  found  that  a  small  number  of  PCA  geometric  modes 
can  capture  manufactured  variability  in  bladed  disks  accurately.29,41,42  Such  reduced  geometric  variability 
models  have  been  used  to  investigate  the  impact  of  geometric  variability  on  axial  compressor  steady  aero¬ 
dynamic  performance  using  Monte  Carlo  simulation  (MCS)  based  on  a  large-scale  nonlinear  CFD  model.29 
Using  MCS  of  a  CFD  model  to  quantify  the  impact  of  geometric  variability  on  unsteady  performance  is  a 
computationally  prohibitive  proposition.  For  example,  if  the  unsteady  analysis  for  one  geometry  takes  one 
minute  to  compute  (a  conservative  estimate),  the  0(50,000)  such  analyses  that  would  be  required  for  a 
MCS  would  take  roughly  one  month  of  CPU  time.  Therefore,  we  desire  to  obtain  a  reduced-order  model 
that  captures  both  unsteady  response  and  variation  over  blade  geometries.  Our  method  combines  the  re¬ 
duced  geometric  variability  model  and  the  adaptive  model  reduction  methodology  of  Algorithm  1  to  obtain 
a  reduced-order  model  that  is  valid  over  a  range  of  forcing  frequencies,  aerodynamic  damping,  and  small 
perturbations  in  blade  geometries,  and  thus  enables  fast  and  accurate  probabilistic  analysis. 

A.  Blade  forced  response  example 

For  the  example  presented  here,  the  flow  is  modeled  using  the  two-dimensional  Euler  equations  written  at  the 
blade  mid-section.  The  average  geometry  of  the  blade  is  shown  in  Figure  1  along  with  the  unstructured  grid 
for  a  single  blade  passage,  which  contains  4292  triangular  elements.  The  Euler  equations  are  discretized  in 
space  with  the  DG  method  described  in  Section  II.  For  the  case  considered  here,  the  incoming  steady-state 
flow  has  a  Mach  number  of  M  =  0.113  and  a  flow  angle  of  (3  =  59°.  Flow  tangency  boundary  conditions 
are  applied  on  the  blade  surfaces.  To  compute  the  steady-state  flow  for  the  nominal  case,  we  exploit  the 
fact  that  the  rotor  is  cyclically  symmetric;  thus,  the  steady  flow  in  each  blade  passage  is  the  same  and  the 
steady-state  solution  can  be  computed  on  a  computational  domain  that  describes  just  a  single  blade  passage. 
Periodic  boundary  conditions  are  applied  on  the  upper  and  lower  boundaries  of  the  grid  to  represent  the 
effects  of  neighboring  blade  passages. 

A  linearized  model  is  derived  for  unsteady  flow  computations  by  assuming  that  the  unsteady  flow  is 
a  small  deviation  from  steady  state  as  described  in  Section  IIA.  An  affine  dependence  of  the  linearized 
system  matrices  on  the  blade  geometries  is  derived  using  the  method  described  in  Section  IIB.  This  leads 
to  a  system  of  the  form  (18),  (19),  where  the  state  vector,  x(f),  contains  the  unknown  perturbation  flow 
quantities  (density,  Cartesian  momentum  components  and  energy).  For  the  DG  formulation,  the  states  are 
the  coefficients  corresponding  to  each  nodal  finite  element  shape  function.  Using  linear  elements,  there  are 
12  degrees  of  freedom  per  element,  giving  a  total  state-space  size  of  n  =  51, 504  states  per  blade  passage.  For 
the  problem  considered  here,  the  forcing  input,  u(f),  describes  the  unsteady  motion  of  each  blade,  which  in 
this  case  is  assumed  to  be  rigid  plunging  motion  (vertical  motion  with  no  rotation).  The  outputs  of  interest, 
y(t),  are  the  unsteady  lift  forces  generated  on  each  blade.  The  initial  perturbation  flow  is  given  by  x°  =  0. 
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Figure  1.  Geometry  and  CFD  mesh  for  a  single  blade  passage. 


B.  Geometric  variability  model 

Geometric  modes  were  computed  using  a  PCA  of  data  modified  from  145  actual  blades,  measured  at  thirteen 
sections  along  the  radial  direction.  The  mid-section  geometries  were  then  extracted.  Thus  the  parameter 
vector  z  contains  the  normally  distributed  random  variables  that  describe  perturbations  in  the  geometry 
of  each  blade  according  to  the  model  (10).  In  Figure  2,  we  consider  a  geometric  model  that  uses  the  two 
dominant  variability  modes,  ns  =  2.  The  figure  shows  the  lift  coefficient,  Cl ,  and  moment  coefficient,  Cm , 
of  a  blade  in  response  to  a  pulse  input  in  plunge  for  a  particular  geometry  that  corresponds  to  z\  =  1.59, 
Z2  =  1.59.  The  response  is  computed  using  the  exact  linearized  CFD  model,  i.e.  the  system  (11),  (12)  and 
the  approximate  linearized  model  (18),  (19)  with  ns  =  2  geometry  modes.  For  reference,  the  response  of  the 
nominal  blade  is  also  shown  in  the  figure.  It  can  be  seen  that  despite  the  small  perturbation  in  geometry,  the 
change  in  lift  and  moment  coefficient  responses  is  significant.  The  approximate  linearized  geometric  model 
captures  the  unsteady  response  accurately. 


Figure  2.  Lift  coefficient,  Cl ,  and  moment  coefficient,  Cm?  in  response  to  a  pulse  input  in  blade  plunge 
displacement  for  the  nominal  geometry  and  a  perturbed  geometry  described  by  two  geometric  PCA  modes 
with  z\  =  1.59,  Z2  =  1.59.  Perturbed  geometry  results  are  computed  with  both  the  exact  and  approximate 
linearized  CFD  model. 
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Table  1  shows  the  error  in  lift  and  moment  outputs  due  to  the  linearized  geometry  approximation  for 
several  different  blade  geometries  with  a  pulse  input  in  plunge.  The  error  e  is  defined  as  the  2-norm  of  the 
difference  between  the  approximate  and  the  exact  linearized  output  as  a  percentage  of  the  change  between 
the  exact  and  the  nominal  output, 


e  = 


\/fof  llye-ya|||cft 

\Jfof  ||ye  -y0Mdt 


x  100%, 


(38) 


where  ye,  ya,  and  yQ  are  respectively  the  exact,  approximate,  and  nominal  outputs.  In  the  table,  ecM 
denotes  the  error  in  moment  coefficient  response,  while  ecL  denotes  the  error  in  lift  coefficient  response. 
These  computations  were  carried  out  over  the  time  horizon  shown  in  Figure  2,  i.e.  with  tf  =  107.  In  general, 
we  expect  the  quality  of  the  approximate  model  to  be  compromised  as  the  size  of  the  geometric  perturbation 
increases.  The  errors  shown  in  Table  1  for  blade  geometries  in  the  tails  of  the  distribution,  i.e.  those  with 
large  geometry  variation,  are  deemed  to  be  acceptable  for  the  probabilistic  application  of  interest  here.  For 
applications  where  greater  accuracy  for  large  geometry  variations  is  important  (for  example,  determining  the 
probability  of  failure  would  require  the  tail  of  the  distribution  to  be  resolved  accurately),  the  results  suggest 
that  the  approximate  linearized  CFD  system  is  not  appropriate.  In  such  cases,  one  might  consider  including 
more  terms  in  (16),  the  Taylor  series  expansion  of  the  CFD  matrices. 


Table  1.  Error  in  approximate  linearized  model  predictions  for  a  pulse  input  in  blade  displacement  for  several 
different  geometries. 


Variability  amplitudes 

ecM  (%) 

ecL(%) 

0!  =  1.59,  22  =  1-59 

5.04 

2.6 

=  1.59,02  =  -1-59 

0.3 

0.1 

0i  =  -1.59,02  =  -1-59 

2.0 

0.8 

0i  =  3.0, 02  =  3.0 

16.6 

9.2 

0i  =  3.0, 02  =  —3.0 

4.1 

2.3 

0i  =  —3.0, 02  =  —3.0 

12.4 

4.7 

C.  Model  reduction 

To  create  a  reduced-order  model  for  use  in  probabilistic  analysis,  the  adaptive  model  reduction  methodology 
of  Algorithm  1  is  applied  to  this  problem.  Results  are  shown  here  for  the  case  of  two  blades  moving  with 
an  interblade  phase  angle  of  180°.  Each  blade  geometry  is  represented  by  two  variability  modes,  giving 
ns  =  4  geometric  parameters  in  this  example.  Applying  the  adaptive  model  reduction  methodology  with 
s  =  10~6  and  with  the  lift  and  moment  coefficients  as  the  outputs  of  interest  yields  a  reduced-order  model 
of  size  nr  =  438  (for  two  blades) .  Algorithm  1  required  five  adaptive  cycles,  over  which  a  total  of  26  Newton 
iterations  were  performed.  Recall  that  in  solving  the  optimization  problem  to  determine  the  next  sample 
point  z*,  the  Newton  step  is  solved  by  CG,  each  iteration  of  which  requires  a  pair  of  state  and  adjoint 
system  solves.  However,  the  coefficient  matrices  of  these  systems  remain  constant  over  the  CG  iterations 
within  a  Newton  step;  therefore,  if  a  direct  solver  is  feasible,  the  coefficient  matrix  needs  to  be  factored  just 
once  for  each  Newton  iteration,  followed  by  triangular  solves  at  each  CG  iteration.  Using  this  strategy,  the 
computational  cost  to  compute  our  reduced  model  was  thus  of  the  order  of  26  full-scale  matrix  factorizations 
(whose  computational  cost  dominates  the  other  calculations).  In  terms  of  CPU  time,  this  corresponds  to  31.4 
hours  on  a  dual  core  personal  computer  with  3.2GHz  Pentium  processor.  We  note  that  for  larger  systems 
(e.g.  three-dimensional  models),  it  would  be  necessary  to  use  iterative  solvers;  in  that  case,  the  cost  would 
scale  with  the  number  of  CG  iterations. 

Although  the  computational  cost  to  perform  the  reduction  is  relatively  high,  we  now  have  a  reduced 
model  of  size  nr  =  438  that  accurately  captures  the  unsteady  response  of  the  original  two-blade  system  with 
n  =  103, 008  states  over  the  range  of  geometries  described  by  the  four  geometric  parameters.  As  an  example 
of  an  application  for  which  this  reduced  model  is  useful,  we  consider  probabilistic  analysis  of  the  system. 
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Specifically,  we  consider  the  impact  of  blade  geometry  variabilities  on  the  work  per  cycle,  which  is  defined  as 
the  integral  of  the  blade  motion  times  the  lift  force  over  one  unsteady  cycle.  A  MCS  was  performed  in  which 
5000  blade  geometries  were  selected  randomly  from  the  given  distributions  for  each  blade.  The  same  5000 
geometries  were  analyzed  using  the  approximate  linearized  CFD  model  and  the  reduced  model.  Figure  3 
shows  the  resulting  probability  density  functions  (PDFs)  of  work  per  cycle  for  the  first  blade,  computed 
using  the  approximate  linearized  CFD  model  and  the  reduced-order  model.  Figure  4  shows  the  PDFs  of 
work  per  cycle  for  the  second  blade.  Table  2  shows  that  the  CPU  time  required  to  compute  the  reduced 
model  MCS  is  a  factor  of  74  times  smaller  than  that  required  for  the  CFD  MCS.  Note  that  the  observed 
speed-up  factor  in  MCS  time  is  less  than  the  relative  reduction  in  the  number  of  states.  This  is  because 
the  CFD  system  matrices  are  sparse  ( A  £  Rraxn  has  2,846,056  nonzero  entries)  while  the  reduced  matrices 
are  not  (Ar  £  has  191^44  nonzero  entries).  Nonetheless,  the  savings  in  computational  time  are 

substantial,  and  more  than  justify  the  offline  cost  required  to  compute  the  reduced  model.  In  practice,  many 
more  than  5000  blade  geometries  are  required  to  obtain  a  converged  MCS;  in  this  case,  the  computational 
cost  of  using  the  CFD  model  becomes  prohibitive.  These  computational  results  were  obtained  on  a  dual 
core  personal  computer  with  3.2GHz  Pentium  processor,  using  a  direct  sparse  solver  for  the  full  model43  and 
Matlab  for  the  reduced  model. 

Table  2  also  compares  the  statistics  of  the  two  distributions.  It  can  be  seen  from  Figure  3,  Figure  4  and 
Table  2  that  the  reduced-order  model  predicts  the  mean,  variance  and  shape  of  the  distribution  of  work 
per  cycle  accurately.  To  further  verify  the  quality  of  the  reduced  model,  we  apply  the  Kolmogorov-Smirnov 
method,44  to  test  whether  the  reduced  work  per  cycle  results  and  the  full  work  per  cycle  results  are  drawn 
from  a  same  distribution.  The  results  show  that  we  cannot  reject  the  hypothesis  that  the  distribution  is  the 
same  at  a  5%  significance  level.  The  probability  that  the  hypothesis  is  true  is  0.9563  for  the  first  blade  and 
0.9999  for  the  second  blade. 


WPC  WPC 


(a)  Full  WPC  for  Blade  1 


(b)  Reduced  WPC  for  Blade  1 


Figure  3.  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of  work  per  cycle  for  Blade  1. 
MCS  results  are  shown  for  5000  blade  geometries.  The  same  geometries  were  analyzed  in  each  case.  Dashed 
line  denotes  the  mean. 

To  further  compare  the  reduced-order  and  CFD  results,  we  pick  four  particular  geometries  corresponding 
to  the  left  tail,  right  tail,  mid-left  and  mid-right  locations  on  the  PDF  of  the  first  blade  as  indicated  by  the 
circles  in  Figure  3(a).  In  Table  3  the  work  per  cycle  is  given  for  these  four  blade  geometries  as  computed 
by  the  exact  CFD  model,  the  approximate  linearized  CFD  model,  and  the  reduced-order  model.  The  table 
shows  that  again  the  approximate  linearized  CFD  is  in  good  agreement  with  the  exact  CFD,  especially  for 
the  mid-left  and  mid-right  cases,  which  have  smaller  variability.  In  addition,  the  effectiveness  of  the  adaptive 
model  reduction  methodology  of  Algorithm  1  can  be  seen  from  the  good  agreement  between  the  approximate 
linearized  CFD  and  the  reduced  results. 


12  of  15 


American  Institute  of  Aeronautics  and  Astronautics 


I  I 


WPC  WPC 

(a)  Full  WPC  for  Blade  2  (b)  Reduced  WPC  for  Blade  2 


Figure  4.  Comparison  of  linearized  CFD  and  reduced-order  model  predictions  of  work  per  cycle  for  Blade  2. 
MCS  results  are  shown  for  5000  blade  geometries.  The  same  geometries  were  analyzed  in  each  case.  Dashed 
line  denotes  the  mean. 


Table  2.  Linearized  CFD  model  and  reduced-order  model  MCS  results.  Work  per  cycle  (WPC)  is  predicted 
for  blade  plunging  motion  at  an  interblade  phase  angle  of  180°  for  5000  randomly  selected  blade  geometries. 


CFD 

Reduced 

Model  size 

103,008 

438 

Computation  time 

347.8  hours 

4.7  hours 

Blade  1  WPC  mean 

-1.8572 

-1.8574 

Blade  1  WPC  variance 

2.7484e-4 

2.778e-4 

Blade  2  WPC  mean 

-1.8581 

-1.8581 

Blade  2  WPC  variance 

2.7887e-4 

2.8044e-4 

Table  3.  Exact  CFD,  approximate  CFD,  and  reduced-order  model  work  per  cycle  prediction  for  the  four 
geometries  indicated  in  Figure  3(a). 


Exact 

Approximate 

Reduced 

Left  tail 

-1.8973 

-1.9056 

-1.9061 

Mid-left 

-1.8637 

-1.8636 

-1.8639 

Mid-right 

-1.8459 

-1.8455 

-1.8457 

Right  tail 

-1.8014 

-1.8086 

-1.8085 
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V.  Conclusions 


The  key  contributions  of  this  paper  are  the  derivation  of  a  linearized  CFD  model  that  permits  the 
effects  of  geometry  variations  to  be  represented  with  an  explicit  affine  function  and  the  development  of  an 
adaptive  sampling  method  to  derive  a  reduced-order  basis  that  spans  both  forcing  input  and  parameter  space. 
Together,  these  contributions  lead  to  a  computationally  tractable  formulation  for  analyzing  the  effects  of 
small  variations  in  geometry  on  unsteady  aerodynamic  response.  The  methodology  was  demonstrated  here 
for  a  problem  that  is  linear  in  the  state  variables  and  affine  in  the  parameter  variables;  however,  the  adaptive 
greedy  sampling  approach  provides  a  general  framework  that  is  applicable  to  nonlinear  problems.  In  the 
general  nonlinear  case,  however,  one  must  address  the  challenge  of  carrying  out  online  reduced-order  model 
computations  in  an  efficient  manner  that  does  not  depend  on  the  large-scale  system  dimension. 
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